\documentclass{article}
\usepackage{lmodern}
\usepackage[a4paper,scale=0.8]{geometry}
\usepackage{listings}
\usepackage[dvipsnames]{xcolor}
\lstset{
    language=C++,
    basicstyle=\ttfamily,
    breaklines=true,
    keywordstyle=\bfseries\color{NavyBlue}, 
    commentstyle=\itshape\color{black!50!white},
    stringstyle=\bfseries\color{PineGreen!90!black},
    columns=flexible,
    numbers=left,
    numbersep=2em,
    numberstyle=\footnotesize,
    frame=single,
    framesep=1em
}

\usepackage{pgf-umlcd}

\usepackage{ctex}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{float}
\usepackage{datetime}

\title{数值分析作业2}
\author{陈乐瑶  3210103710}
\date{\today}

\begin{document}
\maketitle

\part{插值方法的实现}
\begin{tikzpicture}
    \begin{class}[text width=8cm]{Interpolation}{5,0}
        \attribute{+PointNum: int}
        \attribute{+*DiffTable: DifferenceTable}
        \attribute{+inputPoint: vector$<$double$>$}
        \attribute{+inputValue: vector$<$double$>$}
        \operation{+Interpolation(n,x,y)}
        \operation{$\sim$Interpolation()}
        \operation[0]{+get\_DiffTable(): void}
        \operation[0]{+get\_InterpolationValue(x): double}
    \end{class}
    \begin{class}[text width=4cm]{NewtonFormula}{0,-6}
        \inherit{Interpolation}
        \operation{+NewtonFormula(n,x,y)}
        \operation{$\sim$NewtonFormula()}
    \end{class}
    \begin{class}[text width=6cm]{HermiteFormula}{10,-6}
        \inherit{Interpolation}
        \attribute{+inputDerivative: vector$<$double$>$}
        \operation{+HermiteFormula(n,x,y,dy)}
        \operation{$\sim$HermiteFormula()}
    \end{class}
    \begin{class}[text width=5cm]{DifferenceTable}{5,-9}
        \attribute{+head: int}
        \attribute{+tail: int}
        \attribute{+value: double}
        \operation{+DifferenceTable(h,t,v)}
        \operation{$\sim$DifferenceTable()}
        \operation{+get\_head(): int}
        \operation{+get\_tail(): int}
        \operation{+get\_value(): double}
    \end{class}
    \unidirectionalAssociation{NewtonFormula}{}{}{DifferenceTable}
    \unidirectionalAssociation{HermiteFormula}{}{}{DifferenceTable}
\end{tikzpicture}

\begin{itemize}
    \item $DifferenceTable$类用于存储差商表，包含三个成员变量，分别是head，tail和value。如f[x0,x1,x2]，
        则head=0，tail=2，value=f[x0,x1,x2]。
    \item $Interpolation$类中的$get\_DiffTable()$函数用于计算差商表，将差商表对角线元素的值存储在成员
        DiffTable中作为系数。差商表的具体计算规则由讲义给出方法实现，其中两种继承类的实现有所差异。
    \item $Interpolation$类中的$get\_InterpolationValue(x)$函数用于根据给定的x值计算插值结果。
        利用计算差商表后得到的DiffTable中的系数并根据插值公式计算得到。
\end{itemize}


\part{插值方法的应用}
这里我们仅展示部分cpp文件内容，python作图部分较为简单不作展示。
\section{B:Newton Formula}

\begin{lstlisting}
    //define function for the problem
    class function1: public Function{
        public:
            double get_value(const double x)override{return 1.0/(1.0+x*x);};
    }func1;

    for(int n1 = 2; n1 <= 8; n1 += 2){
        //define the input point
        std::vector<double> inputPoint1(n1+1);
        std::vector<double> inputValue1(n1+1);
        for (int i = 0; i <= n1; i++) {
            inputPoint1[i] = -5+10.0*i/n1;
            inputValue1[i] = func1.get_value(inputPoint1[i]);
        }
        //compute the output point
        NewtonFormula newtonSolver(n1+1, inputPoint1, inputValue1);
        for (int i = 0; i <= pointNum; i++) {
            outputPoint1[i] = -5+10.0*i/pointNum;
            outputValue1[i] = newtonSolver.getInterpolationValue(outputPoint1[i]);
        }
        
        //output the result to plot
        string title = "./data/problemB_" + to_string(n1) + ".txt";
        ofstream fout1(title);
        for (int i = 0; i <= pointNum; i++) 
            fout1 << outputPoint1[i] << " " << outputValue1[i] << endl;
        fout1.close();
    } 
\end{lstlisting}

插值结果在python中绘制如下：
\begin{figure}[htbp]
    \includegraphics[width=10cm]{./pic/problemB.png}
    \caption{Newton Formula}
\end{figure}

\section{C:Chebyshev Interpolation}
Chebyshev插值在Newton插值的基础上，对插值点进行了优化，从而减小了插值误差。
\begin{lstlisting}
    //define function for the problem
    class function2: public Function{
        public:
            double get_value(const double x)override{return 1.0/(1.0+25*x*x);};
    }func2;

    for(int n2 = 5; n2 <= 20; n2 += 5){
        //define the input point (Chebyshev nodes)
        std::vector<double> inputPoint2(n2);
        std::vector<double> inputValue2(n2);
        for (int i = 0; i < n2; i++) {
            inputPoint2[i] = cos((2*i+1)*M_PI/(2*n2));//Chebyshev nodes
            inputValue2[i] = func2.get_value(inputPoint2[i]);
        }
        //compute the output point
        NewtonFormula newtonSolver(n2, inputPoint2, inputValue2);
        for (int i = 0; i <= pointNum; i++) {
            outputPoint2[i] = -1+2.0*i/pointNum;
            outputValue2[i] = newtonSolver.getInterpolationValue(outputPoint2[i]);
        }

        //output the result to plot
        string title = "./data/problemC_" + to_string(n2) + ".txt";
        ofstream fout2(title);
        for (int i = 0; i <= pointNum; i++) 
            fout2 << outputPoint2[i] << " " << outputValue2[i] << endl;
        fout2.close();
    }
\end{lstlisting}

插值结果在python中绘制如下：
\begin{figure}[htbp]
    \includegraphics[width=10cm]{./pic/problemC.png}
    \caption{Chebyshev Interpolation}
\end{figure}

\section{D:car problem}
速度作为位移的导数，因此在Hermite插值中，除了输入点和值，还需要输入导数值。

而在计算导数时，采用位移的导数近似值，即$v = \lim_{\epsilon \to 0}\frac{d(t+\epsilon)-d(t-\epsilon)}{2\epsilon}$。

\begin{lstlisting}
    //define the input point
    std::vector<double> inputPoint3={0.0, 3.0, 5.0, 8.0, 13.0};
    std::vector<double> inputValue3_d={0.0, 225.0, 383.0, 623.0, 993.0};
    std::vector<double> inputValue3_v={75.0, 77.0, 80.0, 74.0, 72.0};
    
    //(a): predict position and speed when t = 10
    double eps = 1e-8;
    HermiteFormula hermiteSolver(5, inputPoint3, inputValue3_d, inputValue3_v);
    double outputValue3_10_d = hermiteSolver.getInterpolationValue(10.0);
    double outputValue3_10_v = (hermiteSolver.getInterpolationValue(10.0+eps) 
                          - hermiteSolver.getInterpolationValue(10.0-eps))/(2*eps);
    cout << "when t = 10, d = " << outputValue3_10_d << endl;
    cout << "when t = 10, v = " << outputValue3_10_v << endl;
    
    //(b): whether the velocity exceeds 81
    double maxVelocity = 0, maxVelocityPoint = 0;
    for (int i = 0; i <= pointNum; i++) {
        outputPoint3[i] = i*13.0/pointNum;
        outputValue3_d[i] = hermiteSolver.getInterpolationValue(outputPoint3[i]);
        outputValue3_v[i] = (hermiteSolver.getInterpolationValue(outputPoint3[i]+eps) 
                            - hermiteSolver.getInterpolationValue(outputPoint3[i]-eps))/(2*eps);
        if (outputValue3_v[i] > maxVelocity){
            maxVelocity = outputValue3_v[i];
            maxVelocityPoint = outputPoint3[i];
        }
    }
    if (maxVelocity > 81)
        cout << "The velocity exceeds 81: when t = " << maxVelocityPoint << ", v = " << maxVelocity << endl;
    else
        cout << "The velocity does not exceed 81: when t = " << maxVelocityPoint << ", v = " << maxVelocity << endl;
    
    //output the result to plot
    ofstream fout3("./data/problemD_d.txt");
    ofstream fout4("./data/problemD_v.txt");
    for (int i = 0; i <= pointNum; i++) {
        fout3 << outputPoint3[i] << " " << outputValue3_d[i] << endl;
        fout4 << outputPoint3[i] << " " << outputValue3_v[i] << endl;
    }
    fout3.close();
    fout4.close();
\end{lstlisting}

插值结果在python中绘制如下：
\begin{figure}[htbp]
    \includegraphics[width=10cm]{./pic/problemD.png}
    \caption{car problem}
\end{figure}

\section{E:larvae problem}
用已知两组数据进行Newton插值，分别得到两个样本的平均体重曲线和未来十五天的预测曲线。

\begin{lstlisting}
    //define the input point
    std::vector<double> inputPoint4={0.0, 6.0, 10.0, 13.0, 17.0, 20.0, 28.0};
    std::vector<double> inputValue4_Sp1={6.67, 17.3, 42.7, 37.3, 30.1, 29.3, 28.7};
    std::vector<double> inputValue4_Sp2={6.67, 16.1, 18.9, 15.0, 10.6, 9.44, 8.89};

    //(a): compute the average weight
    std::vector<double> outputPoint4_ave(pointNum+1), outputValue4_ave_Sp1(pointNum+1), outputValue4_ave_Sp2(pointNum+1);
    NewtonFormula newtonSolverE1(7, inputPoint4, inputValue4_Sp1);
    NewtonFormula newtonSolverE2(7, inputPoint4, inputValue4_Sp2);
    for (int i = 0; i <= pointNum; i++) {
        outputPoint4_ave[i] = 28.0*i/pointNum;
        outputValue4_ave_Sp1[i] = newtonSolverE1.getInterpolationValue(outputPoint4_ave[i]);
        outputValue4_ave_Sp2[i] = newtonSolverE2.getInterpolationValue(outputPoint4_ave[i]);
    }

    //output the result to plot
    ofstream fout5("./data/problemE_a1.txt");
    ofstream fout6("./data/problemE_a2.txt");
    for (int i = 0; i <= pointNum; i++) {
        fout5 << outputPoint4_ave[i] << " " << outputValue4_ave_Sp1[i] << endl;
        fout6 << outputPoint4_ave[i] << " " << outputValue4_ave_Sp2[i] << endl;
    }
    fout5.close();
    fout6.close();
    
    //(b): predict whether the samples will die after 15 days
    double outputValue4_day43_Sp1 = newtonSolverE1.getInterpolationValue(28.0+15.0);
    double outputValue4_day43_Sp2 = newtonSolverE2.getInterpolationValue(28.0+15.0);
    cout << "The predicted weight of the first sample after 15 days: " << outputValue4_day43_Sp1 << endl;
    cout << "The predicted weight of the second sample after 15 days: " << outputValue4_day43_Sp2 << endl;
    std::vector<double> outputPoint4_pre(pointNum+1), outputValue4_pre_Sp1(pointNum+1), outputValue4_pre_Sp2(pointNum+1);
    for (int i = 0; i <= pointNum; i++) {
        outputPoint4_pre[i] = (28.0+15.0)*i/pointNum;
        outputValue4_pre_Sp1[i] = newtonSolverE1.getInterpolationValue(outputPoint4_pre[i]);
        outputValue4_pre_Sp2[i] = newtonSolverE2.getInterpolationValue(outputPoint4_pre[i]);
    }

    //output the result to plot
    ofstream fout7("./data/problemE_b1.txt");
    ofstream fout8("./data/problemE_b2.txt");
    for (int i = 0; i <= pointNum; i++) {
        fout7 << outputPoint4_pre[i] << " " << outputValue4_pre_Sp1[i] << endl;
        fout8 << outputPoint4_pre[i] << " " << outputValue4_pre_Sp2[i] << endl;
    }
    fout7.close();
    fout8.close();

\end{lstlisting}

插值结果在python中绘制如下：
\begin{figure}[htbp]
    \begin{minipage}[t]{0.48\textwidth}
        \includegraphics[width=7cm]{./pic/problemE_average.png}
        \caption{average weight}
    \end{minipage}
    \begin{minipage}[t]{0.48\textwidth}
        \includegraphics[width=7cm]{./pic/problemE_predict.png}
        \caption{predict weight}
    \end{minipage}
    \caption{larvae problem}
\end{figure}

\end{document}
